package scu.maqiang.possion;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Possion1DL21 {
	public Possion1DL21(Mesh1L2 mesh) {
		this.mesh = mesh;
		fs = new FES1L21(mesh);
		A = new SRMatrix(fs.GetNdof());
		x = new double[fs.GetNdof()];
		gradient = new double[1][fs.GetNdof()];
		RHS = new double[fs.GetNdof()];
	}
	
	public static void main(String[] args) {
		Mesh1L2 mesh1D = new Mesh1L2();
		mesh1D.line(0, 1, 16);
        //mesh2D.setDomainLabel(regionFunc, null, label);
        
		Possion1DL21 p2dl21 = new Possion1DL21(mesh1D);
		p2dl21.DirichletProblem(BVPType.COMMON, 1.0, 2.0,
				new int[] {1, 2}, 0.0, "MySolutionsPossion1DL2.dat");


		/*
		Mesh1L2Generator mg = new Mesh1L2Generator();
		Mesh1L2 mesh1D = mg.line(2, 8, 60);
		Possion1DL21 p2dl21 = new Possion1DL21(mesh1D);
		ScalarFunc kxx = (xy, label, param) -> 0.1*(5 - 0.6 * xy[0]);
		ScalarFunc source = (xy, label, param) -> 0.03 * (xy[0] - 6) * (xy[0] - 6) * (xy[0] - 6) * (xy[0] - 6);
		ScalarFunc constF = (xy, label, param) -> -1.0;
		p2dl21.GeneralPossionProblem(BVPType.COMMON, kxx, null, source, null,
				                       new int[] {1}, constF, null, 
				                       null, null, null, 
				                       null, null, null, 
				                       null, null, null, "D:\\DCFem1D.dat");

		 */

/*
		Mesh1L2Generator mg = new Mesh1L2Generator();
		Mesh1L2 mesh1D = mg.line(0, 1, 16, 5.0);
		Possion1DL21 p2dl21 = new Possion1DL21(mesh1D);
		ScalarFunc kxx = (xy, label, param) -> 1.0;
		ScalarFunc source = (xy, label, param) -> 10 * Math.exp(-5.0 * xy[0] * xy[0]);
		ScalarFunc fluxLeft = (xy, label, param) -> -1.0;
		ScalarFunc constBd = (xy, label, param) -> 0.0;
		p2dl21.GeneralPossionProblem(BVPType.COMMON, kxx, null, source, null,
				new int[] {2}, constBd, null,
				new int[] {1}, fluxLeft, null,
				null, null, null,
				null, null, null, "D:\\DCFem1D.dat");

 */
		/*
		Mesh1L2Generator mg = new Mesh1L2Generator();
		Mesh1L2 mesh1D = mg.line(0, 1, 32, 2.0);
		Possion1DL21 p2dl21 = new Possion1DL21(mesh1D);
		p2dl21.Helmholtz(1.0, BVPType.COMMON, -87.4, BVPType.CONSIST_MASS,
				new int[]{2}, 0.2,
				new int[]{1}, -1.0,
				"Helmholtz.dat");

		 */
		
	}
	
	
	public void DirichletProblem(BVPType tp, 
			ScalarFunc kxxFunc, double[] paramKxx, 
			ScalarFunc sourceFunc, double[] paramSource, 
			int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			String reportFile) {
		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);

		if (fixedFunc == null) {
			fs.applyBC_MBN(RHS, Direct.All, 0.0, DirichletBoundaryLabel);
		} else {
			fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		}
		solveAndComputeGradientAndExport(reportFile);
	}

	public void DirichletProblem(BVPType tp, double[][] vecKxx, 
			ScalarFunc sourceFunc, double[] paramSource, 
			int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			String reportFile) {
		fs.assembleStiff(vecKxx, tp, A);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);

		if (fixedFunc == null) {
			fs.applyBC_MBN(RHS, Direct.All, 0.0, DirichletBoundaryLabel);
		} else {
			fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		}
		solveAndComputeGradientAndExport(reportFile);
	}

	public void DirichletProblem(BVPType tp, double[][] vecKxx, double[][] vecSource, 
			int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			String reportFile) {
		fs.assembleStiff(vecKxx, tp, A);
		if (vecSource != null) {
			fs.assembleSource(vecSource, tp, RHS);
		}
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);

		if (fixedFunc == null) {
			fs.applyBC_MBN(RHS, Direct.All, 0.0, DirichletBoundaryLabel);
		} else {
			fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		}
		solveAndComputeGradientAndExport(reportFile);
	}
	
	/**
	 * 
	 * @param tp
	 * @param kxx
	 * @param source
	 * @param DirichletBoundaryLabel
	 * @param fixedBV
	 * @param reportFile
	 */
	public void DirichletProblem(BVPType tp, double kxx, double source, int[] DirichletBoundaryLabel,
			double fixedBV, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.All, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(reportFile);
	}

	public void GeneralPossionProblem(BVPType tp, ScalarFunc kxxFunc, double[] paramKxx, ScalarFunc sourceFunc,
			double[] paramSource, int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			int[] NeumannBoundaryLabel, ScalarFunc fluxFunc, double[] paramFlux, int[] RobinBoundaryLabel,
			ScalarFunc robinFluxFunc, double[] paramRobinFlux, 
			ScalarFunc convectFunc, double[] paramConvectFunc, BVPType robinType, String reportFile) {
		
		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
		
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new ScalarFunc[] { fluxFunc }, paramFlux, tp, RHS, NeumannBoundaryLabel);
		}
		if (RobinBoundaryLabel != null) {
			fs.assembleBoundaryMass(new ScalarFunc[] { convectFunc }, paramConvectFunc, robinType, A,
					RobinBoundaryLabel);
			fs.assembleFlux(new ScalarFunc[] { robinFluxFunc }, paramRobinFlux, robinType, RHS, RobinBoundaryLabel);
		}

		if (DirichletBoundaryLabel != null) {
			fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
			if (fixedFunc != null) {
				fs.applyBC_MBN(RHS, Direct.X, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
			} else {
				fs.applyBC_MBN(RHS, Direct.X, 0.0, DirichletBoundaryLabel);
			}
		}
		solveAndComputeGradientAndExport(reportFile);

	}

	public void GeneralPossionProblem(BVPType tp, double kxx, double source, int[] DirichletBoundaryLabel,
			double fixedBV, int[] NeumannBoundaryLabel, double fluxValue, int[] RobinBoundaryLabel,
			double robinFluxValue, double convectValue, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new double[] { fluxValue }, tp, RHS, NeumannBoundaryLabel);
		}
		if (RobinBoundaryLabel != null) {
			fs.assembleBoundaryMass(new double[] { convectValue }, tp, A, RobinBoundaryLabel);
			fs.assembleFlux(new double[] { robinFluxValue }, tp, RHS, RobinBoundaryLabel);
		}

		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(reportFile);
	}
	
	private void solveAndComputeGradientAndExport(String reportFile) {
		//NewIterSSolver solver = new NewIterSSolver(A, false);
		//solver.PCGSSSOR(RHS, x, 1.50, 1);

//		DirectSSolver solver = new DirectSSolver(A);
//		solver.CholeskySolve(RHS, 1);
//		Arrays.fill(x, 0.0);
		DirectSSolver dsolver = new DirectSSolver(A);
		x = dsolver.CholeskySolve(RHS, 1);
		//fs.computeGradient(x, gradient);
		//double[][] gradient = new double[1][];
		gradient[0] = fs.computeGradient(x);
		System.out.println("Solution: ");
		System.out.println("x0 min:" + MVO.min(x) + "\t" + "max: " + MVO.max(x));
		System.out.println("Solution Gradient: ");
		System.out.println("Dx min:" + MVO.min(gradient[0]) + "\t" + "max: " + MVO.max(gradient[0]));
		double[][] mat = new double[][] { x, gradient[0]};
		if (reportFile != null) {
			mesh.toTecplot(reportFile, mat);
		}
	}

	public void Helmholtz(double kxx, BVPType tpStiff, double rho, BVPType tpMass, int[] DirichletBoundaryLabel,
									  double fixedBV, int[] NeumannBoundaryLabel, double fluxValue, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tpStiff, A);
		fs.assembleMass(new double[] { rho }, tpMass, A);
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new double[] { fluxValue }, tpStiff, RHS, NeumannBoundaryLabel);
		}

		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport2(reportFile);
	}

	private void solveAndComputeGradientAndExport2(String reportFile) {
		//NewIterSSolver solver = new NewIterSSolver(A, false);
		//solver.PCGSSSOR(RHS, x, 1.50, 1);

//		DirectSSolver solver = new DirectSSolver(A);
//		solver.CholeskySolve(RHS, 1);
//		Arrays.fill(x, 0.0);
		DirectSSolver dsolver = new DirectSSolver(A);
		x = dsolver.LUSolve(RHS, 1, 1.0e-10);
		fs.computeGradient(x, gradient);
		System.out.println("Solution: ");
		System.out.println("x0 min:" + MVO.min(x) + "\t" + "max: " + MVO.max(x));
		System.out.println("Solution Gradient: ");
		System.out.println("Dx min:" + MVO.min(gradient[0]) + "\t" + "max: " + MVO.max(gradient[0]));
		double[][] mat = new double[][] { x, gradient[0]};
		if (reportFile != null) {
			mesh.toTecplot(reportFile, mat);
		}
	}
	public double[] exportSolution() {
		return x.clone();
	}
	
	public double[][] exportGradient() {
		double[][] gd = new double[1][];
		gd[0] = gradient[0].clone();
		return gd;
	}
	
	Mesh1L2 mesh;
	FES1L21 fs;
	SRMatrix A;
	double[] RHS;
	double[] x;
	double[][] gradient;

}
